* GECM DGP

cscript
* Set working directory
cd "/Users/alikagalwala/Library/CloudStorage/Dropbox/Whitten_Kagalwala/TS/JOP/Acceptance/KW_Replication/GECM_DGP"

do "GECM-DGP.do"

// Monte Carlo Analyses

* Create a matrix to store results
matrix mc_summary_stats = J(7, 72, .)
matrix colnames mc_summary_stats = t phi_u beta_x stat_y  ///
phi_gecm phi_gecm_t1 b1_gecm b1_gecm_t1 b2_gecm b2_gecm_t1 lrm_gecm lrm_gecm_t1 ///
phi_adl phi_adl_t1 b1_adl b1_adl_t1 b2_adl b2_adl_t1 lrm_adl lrm_adl_t1 ///
phi1_adl22 phi1_adl22_t1 phi2_adl22 phi2_adl22_t1 b1_adl22 b1_adl22_t1 b2_adl22 b2_adl22_t1 b3_adl22 b3_adl22_t1 lrm_adl22 lrm_adl22_t1 ///
phi1_adl33 phi1_adl33_t1 phi2_adl33 phi2_adl33_t1 phi3_adl33 phi3_adl33_t1 b1_adl33 b1_adl33_t1 b2_adl33 b2_adl33_t1 b3_adl33 b3_adl33_t1 b4_adl33 b4_adl33_t1 lrm_adl33 lrm_adl33_t1 ///
phi_gecm_power b1_gecm_power b2_gecm_power lrm_gecm_power ///
phi_adl_power b1_adl_power b2_adl_power lrm_adl_power ///
phi1_adl22_power b1_adl22_power b2_adl22_power lrm_adl22_power ///
phi1_adl33_power b1_adl33_power b2_adl33_power lrm_adl33_power ///
ftest_y_adl22_t1 ftest_y_adl22_power ftest_x_adl22_t1 ftest_x_adl22_power ///
ftest_y_adl33_t1 ftest_y_adl33_power ftest_x_adl33_t1 ftest_x_adl33_power
 
* Create a counter to store results in the matrix
local j = 1

foreach t in 25 50 75 100 150 500 1000 {
					parallel sim, expr(stat_y = r(stat_y) phi_gecm = r(phi_gecm) phi_gecm_t1 = r(phi_gecm_t1) ///
					b1_gecm = r(b1_gecm) b1_gecm_t1 = r(b1_gecm_t1) b2_gecm = r(b2_gecm) b2_gecm_t1 = r(b2_gecm_t1) ///
					lrm_gecm = r(lrm_gecm) lrm_gecm_t1 = r(lrm_gecm_t1) ///
					phi_adl = r(phi_adl) phi_adl_t1 = r(phi_adl_t1) b1_adl = r(b1_adl) b1_adl_t1 = r(b1_adl_t1) ///
					b2_adl = r(b2_adl) b2_adl_t1 = r(b2_adl_t1) lrm_adl = r(lrm_adl) lrm_adl_t1 = r(lrm_adl_t1) ///
					phi1_adl22 = r(phi1_adl22) phi1_adl22_t1 = r(phi1_adl22_t1) phi2_adl22 = r(phi2_adl22) phi2_adl22_t1 = r(phi2_adl22_t1) ///
					b1_adl22 = r(b1_adl22) b1_adl22_t1 = r(b1_adl22_t1) b2_adl22 = r(b2_adl22) b2_adl22_t1 = r(b2_adl22_t1) ///
					b3_adl22 = r(b3_adl22) b3_adl22_t1 = r(b3_adl22_t1) lrm_adl22 = r(lrm_adl22) lrm_adl22_t1 = r(lrm_adl22_t1) ///
					phi1_adl33 = r(phi1_adl33) phi1_adl33_t1 = r(phi1_adl33_t1) phi2_adl33 = r(phi2_adl33)  phi2_adl33_t1 = r(phi2_adl33_t1) ///
					phi3_adl33 = r(phi3_adl33) phi3_adl33_t1 = r(phi3_adl33_t1) b1_adl33 = r(b1_adl33) b1_adl33_t1 = r(b1_adl33_t1) b2_adl33 = r(b2_adl33) b2_adl33_t1 = r(b2_adl33_t1) b3_adl33 = r(b3_adl33) b3_adl33_t1 = r(b3_adl33_t1) b4_adl33 = r(b4_adl33) b4_adl33_t1 = r(b4_adl33_t1) ///
					lrm_adl33 = r(lrm_adl33) lrm_adl33_t1 = r(lrm_adl33_t1) ///
					phi_gecm_power = r(phi_gecm_power)  b1_gecm_power = r(b1_gecm_power) b2_gecm_power = r(b2_gecm_power) lrm_gecm_power = r(lrm_gecm_power) ///
					phi_adl_power = r(phi_adl_power) b1_adl_power = r(b1_adl_power) b2_adl_power = r(b2_adl_power) lrm_adl_power = r(lrm_adl_power) ///
					phi1_adl22_power = r(phi1_adl22_power) b1_adl22_power = r(b1_adl22_power) b2_adl22_power = r(b2_adl22_power) lrm_adl22_power = r(lrm_adl22_power) ///
					phi1_adl33_power = r(phi1_adl33_power) b1_adl33_power = r(b1_adl33_power) b2_adl33_power = r(b2_adl33_power) lrm_adl33_power = r(lrm_adl33_power) ftest_y_adl22_t1 = r(ftest_y_adl22_t1) ftest_y_adl22_power = r(ftest_y_adl22_power) ///
					ftest_x_adl22_t1 = r(ftest_x_adl22_t1) ftest_x_adl22_power = r(ftest_x_adl22_power) ftest_y_adl33_t1 = r(ftest_y_adl33_t1) ftest_y_adl33_power = r(ftest_y_adl33_power) ftest_x_adl33_t1 = r(ftest_x_adl33_t1)  ///
					ftest_x_adl33_power = r(ftest_x_adl33_power)) ///
					reps(1000) noisily seeds(12344 12344 12344 12344 12344 12344 12344 12344): ecm, obs(`t') phi_u(0.6) beta_x(1)
			
					matrix mc_summary_stats[`j', 1] = `t'
					matrix mc_summary_stats[`j', 2] = 0.6
					matrix mc_summary_stats[`j', 3] = 1
					
					* Proportion of times y is stationary
					qui su stat_y
					matrix mc_summary_stats[`j', 4] = r(mean)
					
					* GECM model
					qui su phi_gecm
					matrix mc_summary_stats[`j', 5] = r(mean)
					qui su phi_gecm_t1
					matrix mc_summary_stats[`j', 6] = r(mean)
					qui su b1_gecm
					matrix mc_summary_stats[`j', 7] = r(mean)
					qui su b1_gecm_t1
					matrix mc_summary_stats[`j', 8] = r(mean)
					qui su b2_gecm
					matrix mc_summary_stats[`j', 9] = r(mean)
					qui su b2_gecm_t1
					matrix mc_summary_stats[`j', 10] = r(mean)
					qui su lrm_gecm
					matrix mc_summary_stats[`j', 11] = r(mean)
					qui su lrm_gecm_t1
					matrix mc_summary_stats[`j', 12] = r(mean)
			
					* ADL model
					qui su phi_adl
					matrix mc_summary_stats[`j', 13] = r(mean)
					qui su phi_adl_t1
					matrix mc_summary_stats[`j', 14] = r(mean)
					qui su b1_adl
					matrix mc_summary_stats[`j', 15] = r(mean)
					qui su b1_adl_t1
					matrix mc_summary_stats[`j', 16] = r(mean)
					qui su b2_adl
					matrix mc_summary_stats[`j', 17] = r(mean)
					qui su b2_adl_t1
					matrix mc_summary_stats[`j', 18] = r(mean)
					qui su lrm_adl
					matrix mc_summary_stats[`j', 19] = r(mean)
					qui su lrm_adl_t1
					matrix mc_summary_stats[`j', 20] = r(mean)
			
					* ADL (2,2)
					qui su phi1_adl22
					matrix mc_summary_stats[`j', 21] = r(mean)
					qui su phi1_adl22_t1
					matrix mc_summary_stats[`j', 22] = r(mean)
					qui su phi2_adl22
					matrix mc_summary_stats[`j', 23] = r(mean)
					qui su phi2_adl22_t1
					matrix mc_summary_stats[`j', 24] = r(mean)
					qui su b1_adl22
					matrix mc_summary_stats[`j', 25] = r(mean)
					qui su b1_adl22_t1
					matrix mc_summary_stats[`j', 26] = r(mean)
					qui su b2_adl22
					matrix mc_summary_stats[`j', 27] = r(mean)
					qui su b2_adl22_t1
					matrix mc_summary_stats[`j', 28] = r(mean)
					qui su b3_adl22
					matrix mc_summary_stats[`j', 29] = r(mean)
					qui su b3_adl22_t1
					matrix mc_summary_stats[`j', 30] = r(mean)
					qui su lrm_adl22
					matrix mc_summary_stats[`j', 31] = r(mean)
					qui su lrm_adl22_t1
					matrix mc_summary_stats[`j', 32] = r(mean)
					
					
					* ADL(3,3) model
					qui su phi1_adl33
					matrix mc_summary_stats[`j', 33] = r(mean)
					qui su phi1_adl33_t1
					matrix mc_summary_stats[`j', 34] = r(mean)
					qui su phi2_adl33
					matrix mc_summary_stats[`j', 35] = r(mean)
					qui su phi2_adl33_t1
					matrix mc_summary_stats[`j', 36] = r(mean)
					qui su phi3_adl33
					matrix mc_summary_stats[`j', 37] = r(mean)
					qui su phi3_adl33_t1
					matrix mc_summary_stats[`j', 38] = r(mean)
					qui su b1_adl33
					matrix mc_summary_stats[`j', 39] = r(mean)
					qui su b1_adl33_t1
					matrix mc_summary_stats[`j', 40] = r(mean)
					qui su b2_adl33
					matrix mc_summary_stats[`j', 41] = r(mean)
					qui su b2_adl33_t1
					matrix mc_summary_stats[`j', 42] = r(mean)
					qui su b3_adl33
					matrix mc_summary_stats[`j', 43] = r(mean)
					qui su b3_adl33_t1
					matrix mc_summary_stats[`j', 44] = r(mean)
					qui su b4_adl33
					matrix mc_summary_stats[`j', 45] = r(mean)
					qui su b4_adl33_t1
					matrix mc_summary_stats[`j', 46] = r(mean)
					qui su lrm_adl33
					matrix mc_summary_stats[`j', 47] = r(mean)
					qui su lrm_adl33_t1
					matrix mc_summary_stats[`j', 48] = r(mean)
					
					* Power calculations
					qui su phi_gecm_power
					matrix mc_summary_stats[`j', 49] = r(mean)
					qui su b1_gecm_power
					matrix mc_summary_stats[`j', 50] = r(mean)
					qui su b2_gecm_power
					matrix mc_summary_stats[`j', 51] = r(mean)
					qui su lrm_gecm_power
					matrix mc_summary_stats[`j', 52] = r(mean)
					
					qui su phi_adl_power
					matrix mc_summary_stats[`j', 53] = r(mean)
					qui su b1_adl_power
					matrix mc_summary_stats[`j', 54] = r(mean)
					qui su b2_adl_power
					matrix mc_summary_stats[`j', 55] = r(mean)
					qui su lrm_adl_power
					matrix mc_summary_stats[`j', 56] = r(mean)
			
					
					qui su phi1_adl22_power
					matrix mc_summary_stats[`j', 57] = r(mean)
					qui su b1_adl22_power
					matrix mc_summary_stats[`j', 58] = r(mean)
					qui su b2_adl22_power
					matrix mc_summary_stats[`j', 59] = r(mean)
					qui su lrm_adl22_power
					matrix mc_summary_stats[`j', 60] = r(mean)
				
					
					qui su phi1_adl33_power
					matrix mc_summary_stats[`j', 61] = r(mean)
					qui su b1_adl33_power
					matrix mc_summary_stats[`j', 62] = r(mean)
					qui su b2_adl33_power
					matrix mc_summary_stats[`j', 63] = r(mean)
					qui su lrm_adl33_power
					matrix mc_summary_stats[`j', 64] = r(mean)
					
					* F-tests
					qui su ftest_y_adl22_t1
					matrix mc_summary_stats[`j', 65] = r(mean)
					qui su ftest_y_adl22_power
					matrix mc_summary_stats[`j', 66] = r(mean)
					qui su ftest_x_adl22_t1
					matrix mc_summary_stats[`j', 67] = r(mean)
					qui su ftest_x_adl22_power
					matrix mc_summary_stats[`j', 68] = r(mean)
					qui su ftest_y_adl33_t1
					matrix mc_summary_stats[`j', 69] = r(mean)
					qui su ftest_y_adl33_power
					matrix mc_summary_stats[`j', 70] = r(mean)
					qui su ftest_x_adl33_t1
					matrix mc_summary_stats[`j', 71] = r(mean)
					qui su ftest_x_adl33_power
					matrix mc_summary_stats[`j', 72] = r(mean)
					* Iterate the counter
					local j = `j' + 1
		}
		
// Store MC matrix results as a dataframe
clear
set obs 7
svmat double mc_summary_stats, names(col)

scalar gamma1 = 1.5

gen phi_bias_gecm = phi_gecm - (phi_u - 1)
gen b1_bias_gecm = b1_gecm - beta_x
gen b2_bias_gecm = b2_gecm - gamma1
gen lrm_bias_gecm = lrm_gecm - ((-gamma1)/(phi_u -1))

gen phi_bias_adl = phi_adl - phi_u
gen b1_bias_adl = b1_adl - beta_x
gen b2_bias_adl = b2_adl - (gamma1-beta_x)
gen lrm_bias_adl = lrm_adl - ((beta_x + (gamma1-beta_x))/(1-phi_u))

gen phi1_bias_adl22 = phi1_adl22 - phi_u
gen b1_bias_adl22 = b1_adl22 - beta_x
gen b2_bias_adl22 = b2_adl22 - (gamma1-beta_x)
gen lrm_bias_adl22 = lrm_adl22 - ((beta_x + (gamma1-beta_x))/(1-phi_u))
gen phi2_bias_adl22 = phi2_adl22 - 0
gen b3_bias_adl22 = b3_adl22 - 0

gen phi1_bias_adl33 = phi1_adl33 - phi_u
gen phi2_bias_adl33 = phi2_adl33 - 0
gen phi3_bias_adl33 = phi3_adl33 - 0
gen b1_bias_adl33 = b1_adl33 - beta_x
gen b2_bias_adl33 = b2_adl33 - (gamma1-beta_x)
gen b3_bias_adl33 = b3_adl33 - 0
gen b4_bias_adl33 = b4_adl33 - 0
gen lrm_bias_adl33 = lrm_adl33 - ((beta_x + (gamma1-beta_x))/(1-phi_u))

save gecm.dta, replace
